Mosquito microbes analysis
2023-01-31
1 Setup
1.1 Libs
## Loading required package: carData## Loading required package: lattice## Warning: package 'lattice' was built under R version 4.3.1## Loading required package: plyr## Warning: package 'plyr' was built under R version 4.3.1## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 4.3.1##
## Attaching package: 'ggpubr'## The following object is masked from 'package:plyr':
##
## mutate##
## Attaching package: 'cowplot'## The following object is masked from 'package:ggpubr':
##
## get_legend## Loading required package: Matrix## Warning: package 'Matrix' was built under R version 4.3.1## Warning: package 'glmmTMB' was built under R version 4.3.1## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.6.0
## Current Matrix version is 1.6.1.1
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package## Warning: package 'emmeans' was built under R version 4.3.1## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data## Loading required package: MASS##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyser## Loading required package: stats4## This is DHARMa 0.4.6. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')##
## Attaching package: 'nlme'## The following object is masked from 'package:lme4':
##
## lmList##
## Attaching package: 'bestNormalize'## The following object is masked from 'package:MASS':
##
## boxcox1.2 Data
setwd("~/nicolagk@hawaii.edu - Google Drive/My Drive/Mosquito_business/Mosquito_microbes_git/04.microbiome_analysis/04a.diversity")
ps.clean <- readRDS("../../02.process_asvs/ps.clean.trim.rds")
#ps.clean <- readRDS("../../02.process_asvs/ps.clean.trim.rare9200.rds")
samdf <- data.frame(ps.clean@sam_data)
##made below
#df.div <- read.csv("div.trim.csv",row.names=1)2 Alpha diversity calculations (first run through)
##option to run rarefied stuff:
#ps.clean <- ps.rare
df.all <- data.frame(estimate_richness(ps.clean, split=TRUE, measures=c("Shannon","InvSimpson","Observed")))
df.all$sample_name <- rownames(df.all)
df.all.div <- merge(df.all,samdf,by="sample_name") #add sample data
#shannon diversity divided by species richness
df.all.div$even <- df.all.div$Shannon/(log(df.all.div$Observed))
df.div <- df.all.div
#adding most recent sample counts
sums <- data.frame(sample_sums(ps.clean))
colnames(sums) <- c("lib_size_trim")
sums$sample_name <- row.names(sums)
df.div <- merge(df.div,sums,by="sample_name")
#write.csv(df.div,"div.trim.csv")3 Analyzing alpha metrics
3.1 All data
Very strange/interesting pattern where the water types that shouldn’t have a lot of diversity have a lot of OTUs, but not a lot of reads…
ggplot(df.div,aes(x=type,y=totsum,color=infusion))+
geom_boxplot()+
geom_jitter(position=position_jitterdodge())
3.2 Plotz
#just experimental types
df.div.exp <- subset(df.div,type=="Microbial Water"|type=="A.albopictus")
df.div.exp$inf.temp <- paste0(df.div.exp$infusion,df.div.exp$temperature)
df.div.exp$type <- sub("A.albopictus","Mosquitoes",df.div.exp$type)
df.div.exp$type <- sub("Microbial Water","Meso. water",df.div.exp$type)
df.div.exp$type <- factor(df.div.exp$type,levels=c("Mosquitoes","Meso. water"))3.2.1 Richness
df.div.exp.se <- summarySE(df.div.exp,measurevar="Observed",groupvars=c("inf.temp","infusion","temperature","type"))
##dots & error bars
gg.rich <- ggplot(df.div.exp,aes(x=infusion,y=Observed,fill=inf.temp,shape=inf.temp))+
geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
geom_errorbar(data=df.div.exp.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
#geom_boxplot(outlier.shape=NA,alpha=0.5)+
facet_wrap(~type)+
theme_cowplot()+
ylab("ASV richness")+
xlab("Infusion")+
scale_x_discrete(labels=c("OL","SG","PW"))+
scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
scale_shape_manual(values=c(22,24,22,24,22,24),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))
gg.rich
##water by time
df.div.water <- subset(df.div.exp,type=="Meso. water")
gg.mw.rich.time <- ggplot(df.div.water,aes(x=infusion,y=Observed,fill=inf.temp))+
geom_boxplot()+
#geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
#geom_errorbar(data=df.div.exp.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
#geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
#geom_boxplot(outlier.shape=NA,alpha=0.5)+
facet_wrap(~day)+
theme_cowplot()+
ylab("ASV richness")+
xlab("Infusion")+
scale_x_discrete(labels=c("OL","SG","PW"))+
scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
ggtitle("Mesocosm water")3.2.2 Simpson’s
df.div.exp.se <- summarySE(df.div.exp,measurevar="InvSimpson",groupvars=c("inf.temp","infusion","temperature","type"))
##dots & error bars
gg.simp <- ggplot(df.div.exp,aes(x=infusion,y=InvSimpson,fill=inf.temp,shape=inf.temp))+
geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
geom_errorbar(data=df.div.exp.se,aes(ymax=InvSimpson+se,ymin=InvSimpson-se),width=0.2,color="black",position=position_dodge(width=0.6))+
geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
#geom_boxplot(outlier.shape=NA,alpha=0.5)+
facet_wrap(~type)+
theme_cowplot()+
ylab("Simpson's (inv.)")+
xlab("Infusion")+
scale_x_discrete(labels=c("OL","SG","PW"))+
scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
scale_shape_manual(values=c(22,24,22,24,22,24),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))
gg.simp
gg.mw.simp.time <- ggplot(df.div.water,aes(x=infusion,y=InvSimpson,fill=inf.temp))+
geom_boxplot()+
#geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
#geom_errorbar(data=df.div.exp.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
#geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
#geom_boxplot(outlier.shape=NA,alpha=0.5)+
facet_wrap(~day)+
theme_cowplot()+
ylab("Simpson's index (inv.)")+
xlab("Infusion")+
scale_x_discrete(labels=c("OL","SG","PW"))+
scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
ggtitle("Mesocosm water")
3.3 Statz
3.3.1 Richness
3.3.1.1 Water richness stats
df.div.mw <- subset(df.div.exp,type=="Meso. water")
df.div.mw$day <- as.factor(df.div.mw$day)
str(df.div.mw)## 'data.frame': 211 obs. of 32 variables:
## $ sample_name : chr "M_D13_OL1_1_S216_L001" "M_D13_OL1_2_S217_L001" "M_D13_OL1_3_S218_L001" "M_D13_OL1_4_S219_L001" ...
## $ Observed : num 25 28 22 28 31 31 28 29 27 21 ...
## $ Shannon : num 1.63 1.53 1.58 1.53 1.75 ...
## $ InvSimpson : num 2.67 2.4 2.44 2.32 2.93 ...
## $ orgname : chr "M_D13_OL1.1" "M_D13_OL1.2" "M_D13_OL1.3" "M_D13_OL1.4" ...
## $ X16scq : num NA NA NA NA NA NA NA NA NA NA ...
## $ totsum : int 25797 34139 15861 49298 44559 32248 32456 39229 31583 23414 ...
## $ rspsum : int 23466 31716 14492 45548 39491 29021 29261 35577 29041 21719 ...
## $ mesocosm : chr "OL1.1" "OL1.2" "OL1.3" "OL1.4" ...
## $ type : Factor w/ 2 levels "Mosquitoes","Meso. water": 2 2 2 2 2 2 2 2 2 2 ...
## $ stage : chr "Water" "Water" "Water" "Water" ...
## $ sex : chr "" "" "" "" ...
## $ date.collected : chr "8/27/19" "8/27/19" "8/27/19" "8/27/19" ...
## $ day : Factor w/ 3 levels "4","12","20": 2 2 2 2 2 2 2 2 2 2 ...
## $ abdomen.length.mm: num NA NA NA NA NA NA NA NA NA NA ...
## $ wing.length : num NA NA NA NA NA NA NA NA NA NA ...
## $ hindleg.length : num NA NA NA NA NA NA NA NA NA NA ...
## $ special : chr "" "" "" "" ...
## $ dispersal : chr "N" "N" "D" "D" ...
## $ temperature : chr "C" "C" "C" "C" ...
## $ infusion : chr "OL" "OL" "OL" "OL" ...
## $ Wolb_ct1 : num NA NA NA NA NA NA NA NA NA NA ...
## $ Wolb_ct2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ Act_ct : num NA NA NA NA NA NA NA NA NA NA ...
## $ notes : chr "" "" "" "" ...
## $ newday : chr "early" "early" "early" "early" ...
## $ lib_size : num 25797 34139 15861 49298 44559 ...
## $ is.neg : logi FALSE FALSE FALSE FALSE FALSE FALSE ...
## $ lib_size_clean : num 25797 34139 15861 49298 44559 ...
## $ even : num 0.506 0.458 0.511 0.46 0.51 ...
## $ lib_size_trim : num 25797 34134 15861 49290 44524 ...
## $ inf.temp : chr "OLC" "OLC" "OLC" "OLC" ...##
## Shapiro-Wilk normality test
##
## data: df.div.mw$Observed
## W = 0.9775, p-value = 0.001855#wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),family="poisson",data=df.div.mw)
#wrich2<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),family="compois",data=df.div.mw)
wrich3<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
#AICtab(wrich1,wrich2,wrich3)
##full model
wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
##no dispersal
wrich1nd<-glmmTMB(Observed~day+temperature+infusion+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
##minus infusion
wrich1ni<-glmmTMB(Observed~day+temperature+dispersal+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
##minus temperature
wrich1nt<-glmmTMB(Observed~day+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
##minus time
wrich1nnd<-glmmTMB(Observed~temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
Anova(wrich1)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## day 12.9404 2 0.001549 **
## temperature 0.6742 1 0.411577
## infusion 151.4172 2 < 2.2e-16 ***
## dispersal 0.7306 1 0.392704
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Response: Observed
# Chisq Df Pr(>Chisq)
# day 13.2024 2 0.001359 **
# temperature 0.6742 1 0.411599
# infusion 152.9448 2 < 2.2e-16 ***
# dispersal 0.7326 1 0.392033
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
AICtab(wrich1,wrich1nd,wrich1ni,wrich1nt,wrich1nnd) #no temp best... almost tied with no dispersal, then full## dAIC df
## wrich1nt 0.0 8
## wrich1nd 0.1 8
## wrich1 1.3 9
## wrich1nnd 9.9 7
## wrich1ni 82.7 7##interaction?
wrich1.allint <- glmmTMB(Observed~day*temperature*infusion+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.allint1 <- glmmTMB(Observed~day*temperature+infusion+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.allint2 <- glmmTMB(Observed~day+temperature*infusion+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.allint3 <- glmmTMB(Observed~day*infusion+temperature+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw) #second lowest AIC
wrich1.intred <- glmmTMB(Observed~day*infusion+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw) #lowest AIC
wrich1.noint1 <- glmmTMB(Observed~day+infusion+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
wrich1.noint2 <- glmmTMB(Observed~day+infusion+temperature+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
AICtab(wrich1.allint,wrich1.allint1,wrich1.allint2,wrich1.allint3,wrich1.intred,wrich1.noint1,wrich1.noint2)## dAIC df
## wrich1.intred 0.0 11
## wrich1.allint3 2.5 13
## wrich1.allint 5.3 20
## wrich1.noint1 41.8 7
## wrich1.noint2 43.1 8
## wrich1.allint2 43.7 10
## wrich1.allint1 45.5 10## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## day 18.4193 2 0.0001001 ***
## temperature 0.8225 1 0.3644497
## infusion 177.4337 2 < 2.2e-16 ***
## day:temperature 2.4682 2 0.2911031
## day:infusion 62.6715 4 7.957e-13 ***
## temperature:infusion 4.1045 2 0.1284452
## day:temperature:infusion 5.6801 4 0.2243444
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: Observed
# Chisq Df Pr(>Chisq)
# day 18.6821 2 8.775e-05 ***
# temperature 0.8238 1 0.3641
# infusion 180.0249 2 < 2.2e-16 ***
# day:temperature 2.3143 2 0.3144
# day:infusion 62.0786 4 1.060e-12 ***
# temperature:infusion 4.1335 2 0.1266
# day:temperature:infusion 5.5591 4 0.2346
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Anova(wrich1.intred) #all 0.001, no temperature included## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## day 17.390 2 0.0001674 ***
## infusion 164.303 2 < 2.2e-16 ***
## day:infusion 59.023 4 4.653e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## day 17.4811 2 0.00016 ***
## infusion 167.8537 2 < 2.2e-16 ***
## temperature 0.7595 1 0.38348
## dispersal 0.7627 1 0.38248
## day:infusion 58.9948 4 4.717e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#tapply(df.div.mw$log(Observed), df.div.mw$infusion, mean)
#tapply(df.div.mw$log(Observed), df.div.mw$day, mean)
##model checking
shapiro.test(residuals(wrich1.intred)) #awesome##
## Shapiro-Wilk normality test
##
## data: residuals(wrich1.intred)
## W = 0.9921, p-value = 0.3128
#plotResiduals(wrich1.int)
wrich1.intred.resid <- simulateResiduals(fittedModel = wrich1.intred, plot = T)
## Warning in ensurePredictor(simulationOutput, form): DHARMa:::ensurePredictor:
## character string was provided as predictor. DHARMa has converted to factor
## automatically. To remove this warning, please convert to factor before
## attempting to plot with DHARMa.
##posthoc things
##just without dispersal because not plotting it
##also not plotting day but it should be in the model somewhere already
wrich1.allint3.em <- emmeans(wrich1.allint3,~infusion*temperature)
multcomp::cld(wrich1.allint3.em)## infusion temperature emmean SE df lower.CL upper.CL .group
## SW C 19.1 0.351 198 18.5 19.8 1
## SW H 19.5 0.353 198 18.8 20.1 1
## SG C 21.2 0.355 198 20.5 21.9 2
## SG H 21.5 0.357 198 20.8 22.2 2
## OL C 24.7 0.351 198 24.0 25.4 3
## OL H 25.0 0.353 198 24.3 25.7 3
##
## Results are averaged over the levels of: day, dispersal
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.3.3.1.2 Water richness stats - rarefied
Same as above but without offset for library size
df.div.mw <- subset(df.div.exp,type=="Meso. water")
df.div.mw$day <- as.factor(df.div.mw$day)
str(df.div.mw)## 'data.frame': 211 obs. of 32 variables:
## $ sample_name : chr "M_D13_OL1_1_S216_L001" "M_D13_OL1_2_S217_L001" "M_D13_OL1_3_S218_L001" "M_D13_OL1_4_S219_L001" ...
## $ Observed : num 25 28 22 28 31 31 28 29 27 21 ...
## $ Shannon : num 1.63 1.53 1.58 1.53 1.75 ...
## $ InvSimpson : num 2.67 2.4 2.44 2.32 2.93 ...
## $ orgname : chr "M_D13_OL1.1" "M_D13_OL1.2" "M_D13_OL1.3" "M_D13_OL1.4" ...
## $ X16scq : num NA NA NA NA NA NA NA NA NA NA ...
## $ totsum : int 25797 34139 15861 49298 44559 32248 32456 39229 31583 23414 ...
## $ rspsum : int 23466 31716 14492 45548 39491 29021 29261 35577 29041 21719 ...
## $ mesocosm : chr "OL1.1" "OL1.2" "OL1.3" "OL1.4" ...
## $ type : Factor w/ 2 levels "Mosquitoes","Meso. water": 2 2 2 2 2 2 2 2 2 2 ...
## $ stage : chr "Water" "Water" "Water" "Water" ...
## $ sex : chr "" "" "" "" ...
## $ date.collected : chr "8/27/19" "8/27/19" "8/27/19" "8/27/19" ...
## $ day : Factor w/ 3 levels "4","12","20": 2 2 2 2 2 2 2 2 2 2 ...
## $ abdomen.length.mm: num NA NA NA NA NA NA NA NA NA NA ...
## $ wing.length : num NA NA NA NA NA NA NA NA NA NA ...
## $ hindleg.length : num NA NA NA NA NA NA NA NA NA NA ...
## $ special : chr "" "" "" "" ...
## $ dispersal : chr "N" "N" "D" "D" ...
## $ temperature : chr "C" "C" "C" "C" ...
## $ infusion : chr "OL" "OL" "OL" "OL" ...
## $ Wolb_ct1 : num NA NA NA NA NA NA NA NA NA NA ...
## $ Wolb_ct2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ Act_ct : num NA NA NA NA NA NA NA NA NA NA ...
## $ notes : chr "" "" "" "" ...
## $ newday : chr "early" "early" "early" "early" ...
## $ lib_size : num 25797 34139 15861 49298 44559 ...
## $ is.neg : logi FALSE FALSE FALSE FALSE FALSE FALSE ...
## $ lib_size_clean : num 25797 34139 15861 49298 44559 ...
## $ even : num 0.506 0.458 0.511 0.46 0.51 ...
## $ lib_size_trim : num 25797 34134 15861 49290 44524 ...
## $ inf.temp : chr "OLC" "OLC" "OLC" "OLC" ...##
## Shapiro-Wilk normality test
##
## data: df.div.mw$Observed
## W = 0.9775, p-value = 0.001855#wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),family="poisson",data=df.div.mw)
#wrich2<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),family="compois",data=df.div.mw)
wrich3<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
#AICtab(wrich1,wrich2,wrich3)
##full model
wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
##no dispersal
wrich1nd<-glmmTMB(Observed~day+temperature+infusion+(1|mesocosm), data=df.div.mw)
##minus infusion
wrich1ni<-glmmTMB(Observed~day+temperature+dispersal+(1|mesocosm), data=df.div.mw)
##minus temperature
wrich1nt<-glmmTMB(Observed~day+infusion+dispersal+(1|mesocosm), data=df.div.mw)
##minus time
wrich1nnd<-glmmTMB(Observed~temperature+infusion+dispersal+(1|mesocosm), data=df.div.mw)
Anova(wrich1)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## day 9.4355 2 0.008935 **
## temperature 0.4517 1 0.501514
## infusion 124.7413 2 < 2.2e-16 ***
## dispersal 0.6137 1 0.433416
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: Observed
# Chisq Df Pr(>Chisq)
# day 10.2547 2 0.005932 **
# temperature 0.5020 1 0.478604
# infusion 133.7700 2 < 2.2e-16 ***
# dispersal 0.5322 1 0.465695
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##same p-values as not rarefied
AICtab(wrich1,wrich1nd,wrich1ni,wrich1nt,wrich1nnd) #no temp best... tied with no dispersal, then full## dAIC df
## wrich1nt 0.0 8
## wrich1nd 0.2 8
## wrich1 1.5 9
## wrich1nnd 6.8 7
## wrich1ni 74.9 7##interaction?
wrich1.allint <- glmmTMB(Observed~day*temperature*infusion+(1|mesocosm),data=df.div.mw)
wrich1.allint1 <- glmmTMB(Observed~day*temperature+infusion+(1|mesocosm),data=df.div.mw)
wrich1.allint2 <- glmmTMB(Observed~day+temperature*infusion+(1|mesocosm),data=df.div.mw)
wrich1.allint3 <- glmmTMB(Observed~day*infusion+temperature+dispersal+(1|mesocosm),data=df.div.mw) #second lowest AIC
wrich1.intred <- glmmTMB(Observed~day*infusion+(1|mesocosm), data=df.div.mw) #lowest AIC
wrich1.noint1 <- glmmTMB(Observed~day+infusion+(1|mesocosm), data=df.div.mw)
wrich1.noint2 <- glmmTMB(Observed~day+infusion+temperature+(1|mesocosm), data=df.div.mw)
AICtab(wrich1.allint,wrich1.allint1,wrich1.allint2,wrich1.allint3,wrich1.intred,wrich1.noint1,wrich1.noint2) #intred best, same as not rarfied## dAIC df
## wrich1.intred 0.0 11
## wrich1.allint3 2.8 13
## wrich1.allint 5.0 20
## wrich1.noint1 45.5 7
## wrich1.noint2 47.0 8
## wrich1.allint2 47.4 10
## wrich1.allint1 49.8 10## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## day 13.4995 2 0.001171 **
## temperature 0.5746 1 0.448450
## infusion 153.0372 2 < 2.2e-16 ***
## day:temperature 1.9400 2 0.379090
## day:infusion 67.7062 4 6.919e-14 ***
## temperature:infusion 4.5128 2 0.104729
## day:temperature:infusion 6.4005 4 0.171168
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: Observed
# Chisq Df Pr(>Chisq)
# day 14.7390 2 0.0006302 ***
# temperature 0.6319 1 0.4266476
# infusion 164.6677 2 < 2.2e-16 ***
# day:temperature 2.1143 2 0.3474442
# day:infusion 68.9507 4 3.78e-14 ***
# temperature:infusion 4.6235 2 0.0990855 .
# day:temperature:infusion 6.1549 4 0.1878773
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Anova(wrich1.intred) #all 0.001, no temperature included## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## day 12.743 2 0.001709 **
## infusion 141.367 2 < 2.2e-16 ***
## day:infusion 63.747 4 4.725e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## day 12.8152 2 0.001649 **
## infusion 143.7560 2 < 2.2e-16 ***
## temperature 0.5250 1 0.468715
## dispersal 0.6716 1 0.412511
## day:infusion 63.7223 4 4.781e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#tapply(df.div.mw$log(Observed), df.div.mw$infusion, mean)
#tapply(df.div.mw$log(Observed), df.div.mw$day, mean)
##model checking
shapiro.test(residuals(wrich1.intred)) #awesome##
## Shapiro-Wilk normality test
##
## data: residuals(wrich1.intred)
## W = 0.99008, p-value = 0.1556
#plotResiduals(wrich1.int)
wrich1.intred.resid <- simulateResiduals(fittedModel = wrich1.intred, plot = T)
## Warning in ensurePredictor(simulationOutput, form): DHARMa:::ensurePredictor:
## character string was provided as predictor. DHARMa has converted to factor
## automatically. To remove this warning, please convert to factor before
## attempting to plot with DHARMa.
##posthoc things
##just without dispersal because not plotting it
##also not plotting day but it should be in the model somewhere already
wrich1.allint3.em <- emmeans(wrich1.allint3,~infusion*temperature)
multcomp::cld(wrich1.allint3.em)## infusion temperature emmean SE df lower.CL upper.CL .group
## SW C 19.2 0.380 198 18.4 19.9 1
## SW H 19.5 0.382 198 18.7 20.2 12
## SG C 21.0 0.383 198 20.2 21.7 23
## SG H 21.2 0.386 198 20.5 22.0 3
## OL C 24.7 0.380 198 23.9 25.4 4
## OL H 25.0 0.382 198 24.2 25.7 4
##
## Results are averaged over the levels of: day, dispersal
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.3.3.1.3 Mosquito richness stats
##
## Shapiro-Wilk normality test
##
## data: df.div.mq$Observed
## W = 0.9207, p-value = 9.183e-09##different stat families
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
mrich2 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="compois",data=df.div.mq)
mrich3 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="genpois",data=df.div.mq)## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluationmrich4 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="poisson",data=df.div.mq)
#mrich4 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="nbinom1",data=df.div.mq)
#mrich5 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="nbinom2",data=df.div.mq)
AICtab(mrich1,mrich2,mrich3,mrich4) #gaussian again## dAIC df
## mrich1 0.0 10
## mrich2 62.3 10
## mrich3 63.9 10
## mrich4 64.6 9##full model
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##no dispersal
mrich1nd <- glmmTMB(Observed~newday+temperature+infusion+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##minus infusion
mrich1ni <- glmmTMB(Observed~newday+temperature+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##minus temperature
mrich1nt <- glmmTMB(Observed~newday+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##minus time
mrich1nnd <- glmmTMB(Observed~temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mq)
##minus sex
mrich1ns <- glmmTMB(Observed~newday+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
Anova(mrich1)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## newday 2.8486 2 0.2407
## temperature 0.0704 1 0.7908
## infusion 0.6272 2 0.7308
## dispersal 0.9250 1 0.3362
## sex 0.0527 1 0.8185##
## Shapiro-Wilk normality test
##
## data: residuals(mrich1)
## W = 0.93556, p-value = 1.302e-07
## infusion temperature emmean SE df lower.CL upper.CL .group
## OL H 8.05 0.395 185 7.27 8.83 1
## OL C 8.17 0.323 185 7.53 8.80 1
## SG H 8.19 0.425 185 7.35 9.03 1
## SG C 8.31 0.465 185 7.39 9.23 1
## SW H 8.51 0.454 185 7.61 9.40 1
## SW C 8.63 0.531 185 7.58 9.67 1
##
## Results are averaged over the levels of: newday, dispersal, sex
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.3.3.1.4 Mosquito richness stats - rarefied
Same as above but without offset for library size
##
## Shapiro-Wilk normality test
##
## data: df.div.mq$Observed
## W = 0.9207, p-value = 9.183e-09##different stat families
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),data=df.div.mq)
mrich2 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="compois",data=df.div.mq)
mrich3 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="genpois",data=df.div.mq)
mrich4 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="poisson",data=df.div.mq)
#mrich4 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="nbinom1",data=df.div.mq)
#mrich5 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="nbinom2",data=df.div.mq)
AICtab(mrich1,mrich2,mrich3,mrich4) #compois ## dAIC df
## mrich2 0.0 10
## mrich3 6.5 10
## mrich4 18.3 9
## mrich1 18.9 10##full model
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="compois",data=df.div.mq)
##no dispersal
mrich1nd <- glmmTMB(Observed~newday+temperature+infusion+sex+(1|mesocosm),family="compois",data=df.div.mq)
##minus infusion
mrich1ni <- glmmTMB(Observed~newday+temperature+dispersal+sex+(1|mesocosm),family="compois",data=df.div.mq)
##minus temperature
mrich1nt <- glmmTMB(Observed~newday+infusion+dispersal+sex+(1|mesocosm),family="compois",data=df.div.mq)
##minus time
mrich1nnd <- glmmTMB(Observed~temperature+infusion+dispersal+sex+(1|mesocosm),family="compois", data=df.div.mq)
##minus sex
mrich1ns <- glmmTMB(Observed~newday+temperature+infusion+dispersal+(1|mesocosm),family="compois",data=df.div.mq)
Anova(mrich1)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## newday 3.9718 2 0.1373
## temperature 0.0016 1 0.9677
## infusion 0.2076 2 0.9014
## dispersal 0.6153 1 0.4328
## sex 0.0050 1 0.9435##
## Shapiro-Wilk normality test
##
## data: residuals(mrich1)
## W = 0.93647, p-value = 1.55e-07
## infusion temperature emmean SE df asymp.LCL asymp.UCL .group
## OL C 2.09 0.0428 Inf 2.00 2.17 1
## OL H 2.09 0.0504 Inf 1.99 2.19 1
## SG C 2.11 0.0606 Inf 1.99 2.23 1
## SG H 2.11 0.0558 Inf 2.00 2.22 1
## SW C 2.11 0.0653 Inf 1.99 2.24 1
## SW H 2.12 0.0573 Inf 2.00 2.23 1
##
## Results are averaged over the levels of: newday, dispersal, sex
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.3.3.2 Simpson’s
3.3.2.1 Water Simpson’s stats
##
## Shapiro-Wilk normality test
##
## data: df.div.mw$InvSimpson
## W = 0.96756, p-value = 9.003e-05
##
## Shapiro-Wilk normality test
##
## data: log(df.div.mw$InvSimpson)
## W = 0.95673, p-value = 5.166e-06## Best Normalizing transformation with 211 Observations
## Estimated Normality Statistics (Pearson P / df, lower => more normal):
## - arcsinh(x): 1.5829
## - Box-Cox: 1.3786
## - Center+scale: 1.3483
## - Double Reversed Log_b(x+a): 1.498
## - Exp(x): 3.4691
## - Log_b(x+a): 1.6639
## - orderNorm (ORQ): 1.1158
## - sqrt(x + a): 1.3829
## - Yeo-Johnson: 1.385
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##
## Based off these, bestNormalize chose:
## orderNorm Transformation with 211 nonmissing obs and no ties
## - Original quantiles:
## 0% 25% 50% 75% 100%
## 1.359 1.947 2.669 3.225 5.130##
## Shapiro-Wilk normality test
##
## data: df.div.mw$simp.norm
## W = 0.96756, p-value = 9.003e-05##just replaced InvSimpson with simp.norm
##full model
wsimp1 <- glmmTMB(simp.norm~day+temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
##no dispersal
wsimp1nd <- glmmTMB(simp.norm~day+temperature+infusion+(1|mesocosm),data=df.div.mw)
##minus infusion
wsimp1ni <- glmmTMB(simp.norm~day+temperature+dispersal+(1|mesocosm),data=df.div.mw)
##minus temperature
wsimp1nt <- glmmTMB(simp.norm~day+infusion+dispersal+(1|mesocosm),data=df.div.mw)
##minus time
wsimp1nnd <- glmmTMB(simp.norm~temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
Anova(wsimp1)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: simp.norm
## Chisq Df Pr(>Chisq)
## day 45.4454 2 1.354e-10 ***
## temperature 48.9643 1 2.607e-12 ***
## infusion 276.6117 2 < 2.2e-16 ***
## dispersal 0.0291 1 0.8644
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: simp.norm
# Chisq Df Pr(>Chisq)
# day 49.7901 2 1.542e-11 ***
# temperature 41.3001 1 1.306e-10 ***
# infusion 229.1055 2 < 2.2e-16 ***
# dispersal 0.0091 1 0.924
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##~same results with normalized results
##same results with rarefied
AICtab(wsimp1,wsimp1nd,wsimp1ni,wsimp1nt,wsimp1nnd)## dAIC df
## wsimp1nd 0.0 8
## wsimp1 2.0 9
## wsimp1nt 37.3 8
## wsimp1nnd 37.5 7
## wsimp1ni 112.4 7##no dispersal good, then full model
##interaction?
wsimp1.int1 <- glmmTMB(simp.norm~day*infusion*temperature+(1|mesocosm), data=df.div.mw)
wsimp1.int2 <- glmmTMB(simp.norm~day*infusion+temperature+(1|mesocosm), data=df.div.mw) #best AIC
wsimp1.int3 <- glmmTMB(simp.norm~day+infusion*temperature+(1|mesocosm), data=df.div.mw)
AICtab(wsimp1.int1,wsimp1.int2,wsimp1.int3,wsimp1nd) #int better## dAIC df
## wsimp1.int2 0.0 12
## wsimp1.int1 2.1 20
## wsimp1.int3 33.6 10
## wsimp1nd 36.8 8## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: simp.norm
## Chisq Df Pr(>Chisq)
## day 62.772 2 2.340e-14 ***
## infusion 274.807 2 < 2.2e-16 ***
## temperature 47.216 1 6.358e-12 ***
## day:infusion 53.210 4 7.703e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Response: simp.norm
# Chisq Df Pr(>Chisq)
# day 70.733 2 4.370e-16 ***
# infusion 227.491 2 < 2.2e-16 ***
# temperature 39.606 1 3.107e-10 ***
# day:infusion 58.541 4 5.874e-12 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##all of the below was bad before normalizing
shapiro.test(residuals(wsimp1))##
## Shapiro-Wilk normality test
##
## data: residuals(wsimp1)
## W = 0.95063, p-value = 1.211e-06
## infusion temperature emmean SE df lower.CL upper.CL .group
## SG C -1.298 0.0882 202 -1.472 -1.1239 1
## SG H -0.682 0.0887 202 -0.857 -0.5074 2
## OL C -0.109 0.0875 202 -0.282 0.0635 3
## SW C 0.464 0.0875 202 0.292 0.6368 4
## OL H 0.507 0.0880 202 0.333 0.6801 4
## SW H 1.080 0.0880 202 0.906 1.2533 5
##
## Results are averaged over the levels of: day, dispersal
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same. # infusion temperature emmean SE df lower.CL upper.CL .group
# SG C -1.255 0.0934 202 -1.439 -1.0708 1
# SG H -0.656 0.0939 202 -0.841 -0.4712 2
# OL C -0.104 0.0927 202 -0.287 0.0786 3
# SW C 0.442 0.0927 202 0.259 0.6244 4
# OL H 0.495 0.0932 202 0.311 0.6783 4
# SW H 1.040 0.0932 202 0.857 1.2240 53.3.2.2 Mosquito Simpson’s
From manuscript pre-revisions: “In contrast, the Simpson’s index, which integrates evenness into a measure of alpha diversity, of the adult mosquito microbiome was lower in male mosquitoes relative to female mosquitoes (p<0.0001). This trend was associated with the increased dominance of wAlbB in males. The Simpson’s index of the adult mosquito microbiome did not vary significantly with the dispersal, aquatic chemistry, temperature, or time of emergence.”
##full model
msimp.all <- glmmTMB(simp.norm~newday+temperature+infusion+dispersal+sex+(1|mesocosm), data=df.div.mq)
##without sex
msimp.nosex <- glmmTMB(simp.norm~day+temperature+infusion+dispersal+(1|mesocosm), data=df.div.mq)
##without dispersal
msimp.nodis <- glmmTMB(simp.norm~newday+temperature+infusion+sex+(1|mesocosm), data=df.div.mq)
##without infusion
msimp.noinf <- glmmTMB(simp.norm~newday+temperature+dispersal+sex+(1|mesocosm), data=df.div.mq)
##without temperature
msimp.notem <- glmmTMB(simp.norm~newday+infusion+dispersal+sex+(1|mesocosm), data=df.div.mq)
##without time
msimp.notim <- glmmTMB(simp.norm~temperature+infusion+dispersal+sex+(1|mesocosm), data=df.div.mq)
Anova(msimp.all)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: simp.norm
## Chisq Df Pr(>Chisq)
## newday 2.3243 2 0.31281
## temperature 0.2966 1 0.58605
## infusion 2.7471 2 0.25321
## dispersal 2.8196 1 0.09312 .
## sex 67.6024 1 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: simp.norm
# Chisq Df Pr(>Chisq)
# newday 2.2911 2 0.31805
# temperature 0.2987 1 0.58469
# infusion 2.7011 2 0.25909
# dispersal 2.8332 1 0.09234 .
# sex 67.5319 1 < 2e-16 ***
AICtab(msimp.all,msimp.nosex,msimp.nodis,msimp.noinf,msimp.notem,msimp.notim)## dAIC df
## msimp.notem 0.0 9
## msimp.notim 0.0 8
## msimp.noinf 0.4 8
## msimp.all 1.7 10
## msimp.nodis 2.4 9
## msimp.nosex 59.1 8#no time & no temp tied for best, then no infusion
##all of the below was bad before normalizing
shapiro.test(residuals(msimp.all))##
## Shapiro-Wilk normality test
##
## data: residuals(msimp.all)
## W = 0.98795, p-value = 0.09731
##posthoc things
msimp.all.em <- emmeans(msimp.all,~infusion*temperature)
multcomp::cld(msimp.all.em)## infusion temperature emmean SE df lower.CL upper.CL .group
## SG H -0.2061 0.152 185 -0.506 0.0939 1
## SG C -0.1164 0.167 185 -0.446 0.2136 1
## OL H -0.0764 0.143 185 -0.359 0.2063 1
## OL C 0.0132 0.115 185 -0.215 0.2409 1
## SW H 0.1483 0.162 185 -0.172 0.4682 1
## SW C 0.2379 0.193 185 -0.142 0.6182 1
##
## Results are averaged over the levels of: newday, dispersal, sex
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.3.4 Infusion water things
df.div.inf <- subset(df.div,type=="Infusion water")
df.div.inf.se <- summarySE(df.div.inf,measurevar="Observed",groupvars=c("infusion"))
df.div.inf.se## infusion N Observed sd se ci
## 1 OL 3 44.666667 4.1633320 2.4037009 10.342290
## 2 SG 3 20.666667 2.0816660 1.2018504 5.171145
## 3 SW 3 3.333333 0.5773503 0.3333333 1.434218# gg.iw.rich <- ggplot(df.div.inf,aes(x=infusion,y=Observed,fill=infusion))+
# geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
# geom_errorbar(data=df.div.inf.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
# geom_point(data=df.div.inf.se,size=2.5,position=position_dodge(width=0.6))+
# #geom_boxplot(outlier.shape=NA,alpha=0.5)+
# theme_cowplot()+
# ggtitle("Infusion water")
# gg.iw.rich4 Diversity correlations
4.1 Diversity x day
Doesn’t look interesting
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
4.2 Diversity - mesocosm vs. mosquito
##
## Attaching package: 'reshape'## The following object is masked from 'package:Matrix':
##
## expand## The following object is masked from 'package:cowplot':
##
## stamp## The following objects are masked from 'package:plyr':
##
## rename, round_anymeso.corr <- merge(df.div.mq,df.div.mw,by="mesocosm")
ggplot(meso.corr,aes(x=InvSimpson.x,y=InvSimpson.y))+
geom_point()+
#facet_wrap(~day.y)+
geom_smooth()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'